Chapter 8: Multi-parameter evidence synthesis
Bob
2025-02-28
library('rjags')
library('coda')
library('jagshelper')
library('jagsUI')
library('mcmcplots') # caterplot
library('dplyr')maple_syrup_example
CONFIDENCE PROFILE EXAMPLE: Maple Syrup Urine Disease
Objective: To estimate the probability of retardation for a case of MSUD without screening (theta.n), & change in retardation rate associated with screening (e.d)
# Model
model_code <- "model
{
#All data assumed to arise from binomial distributions with the appropriate parameters
r.r ~ dbin(r, n.r)
r.s ~ dbin(phi.s, n.s)
r.n ~ dbin(phi.n, n.n)
r.em ~ dbin(theta.em, n.em)
r.lm ~ dbin(theta.lm, n.lm)
#Define functional relationships
# Prob retardation for a case of MSUD who is screened
theta.sm <- phi.s * theta.em + (1 - phi.s)*theta.lm
# Prob retardation for a case of MSUD who is NOT screened
theta.nm <- phi.n * theta.em + (1 - phi.n) * theta.lm
# Expected retardation per 100000 newborns who are screened
theta.s <- (theta.sm * r) * 100000
# Expected retardation per 100000 newborns who are NOT screened
theta.n <- (theta.nm * r) * 100000
# Change in retardation rate associated with screening
e.d <- theta.s - theta.n
# Prior distributions - 'non-informative' Beta(1,1) priors
r ~ dbeta( 1, 1)
phi.s ~ dbeta( 1, 1)
phi.n ~ dbeta( 1, 1)
theta.em ~ dbeta( 1, 1)
theta.lm ~ dbeta( 1, 1)
}
" %>% textConnection
# Data
data <- list(
n.r = 724262,
r.r = 7,
n.s = 276,
r.s = 253,
n.n = 18,
r.n = 8,
n.em=10,
r.em = 2,
n.lm = 10,
r.lm = 10)
# Starting/initial values
parameters_to_save <- c("e.d") # What else?
results <- jags(data = data,
# inits = initial_values,
parameters.to.save = parameters_to_save,
model.file = model_code,
n.chains = 1, # length(initial_values),
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000,
n.thin = 2)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 5
## Unobserved stochastic nodes: 5
## Total graph size: 30
##
## Initializing model
##
## Adaptive phase, 100 iterations x 1 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
##
## Calculating statistics.......
##
## Done.summary(results)## mean sd 2.5% 25% 50% 75%
## e.d -0.3411557 0.1665444 -0.7418316 -0.4318169 -0.3127903 -0.2207268
## deviance 20.1054097 3.2767312 15.6809410 17.6874476 19.4715124 21.8245152
## 97.5% Rhat n.eff overlap0 f
## e.d -0.10076 NA 1 0 0.9999333
## deviance 28.26577 NA 1 0 1.0000000plot(results)
hivepi_6
Simple version of the HIV Epi model with just 6 data points.
# Model
model_code <- "model{
# SET PRIORS
a ~ dbeta( 1,2)
z ~ dbeta (1,1)
b <- z * (1-a) # sets constraint (1-a-b > 0)
c ~ dbeta (1,1)
d ~ dbeta (1,1)
e ~ dbeta (1,1)
# VECTOR p[1:6] HOLDS THE EXPECTED PROBABILITIES FOR EACH DATA POINT
p[1] <- a
p[2] <- b
p[3] <- c
p[4] <- d
p[5] <- (d*b + e*(1-a-b))/(1- a)
p[6] <- c*a + d*b + e*(1-a-b)
# LIKELIHOOD AND DIAGNOSTICS
for(i in 1: 6) {
r[i] ~ dbin(p[i],n[i])
rhat[i] <- p[i] * n[i]
dev[i] <- 2 * (r[i] * log(r[i]/rhat[i]) + (n[i]-r[i]) * log((n[i]-r[i])/(n[i]-rhat[i])))
}
resdev <- sum(dev[])
}
" %>% textConnection
# Data
data <- list(
r=c(11044, 12, 252, 10, 74, 254),
n=c(104577, 882, 15428, 473, 136139, 102287)
)
# Starting/initial values: None specified
parameters_to_save <- c("a", "b", "c", "d", "e", "p", "r", "rhat", "dev", "resdev")
results <- jags(data = data,
# inits = initial_values,
parameters.to.save = parameters_to_save,
model.file = model_code,
n.chains = 1, # How many chains?
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000,
n.thin = 2)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 6
## Unobserved stochastic nodes: 5
## Total graph size: 96
##
## Initializing model
##
## Adaptive phase, 100 iterations x 1 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
##
## Calculating statistics.......
##
## Done.summary(results)## mean sd 2.5% 25% 50%
## a 1.057108e-01 9.478321e-04 1.038570e-01 1.050640e-01 1.057094e-01
## b 1.404450e-02 3.680695e-03 7.651316e-03 1.142074e-02 1.375050e-02
## c 1.715934e-02 8.759215e-04 1.547653e-02 1.656171e-02 1.714567e-02
## d 2.186721e-02 6.112196e-03 1.140236e-02 1.751401e-02 2.137004e-02
## e 2.456504e-04 1.226242e-04 2.264500e-05 1.533394e-04 2.460119e-04
## p[1] 1.057108e-01 9.478321e-04 1.038570e-01 1.050640e-01 1.057094e-01
## p[2] 1.404450e-02 3.680695e-03 7.651316e-03 1.142074e-02 1.375050e-02
## p[3] 1.715934e-02 8.759215e-04 1.547653e-02 1.656171e-02 1.714567e-02
## p[4] 2.186721e-02 6.112196e-03 1.140236e-02 1.751401e-02 2.137004e-02
## p[5] 5.799799e-04 6.292478e-05 4.629310e-04 5.360930e-04 5.776298e-04
## p[6] 2.332568e-03 9.657936e-05 2.146206e-03 2.265966e-03 2.331312e-03
## r[1] 1.104400e+04 0.000000e+00 1.104400e+04 1.104400e+04 1.104400e+04
## r[2] 1.200000e+01 0.000000e+00 1.200000e+01 1.200000e+01 1.200000e+01
## r[3] 2.520000e+02 0.000000e+00 2.520000e+02 2.520000e+02 2.520000e+02
## r[4] 1.000000e+01 0.000000e+00 1.000000e+01 1.000000e+01 1.000000e+01
## r[5] 7.400000e+01 0.000000e+00 7.400000e+01 7.400000e+01 7.400000e+01
## r[6] 2.540000e+02 0.000000e+00 2.540000e+02 2.540000e+02 2.540000e+02
## rhat[1] 1.105492e+04 9.912143e+01 1.086105e+04 1.098728e+04 1.105477e+04
## rhat[2] 1.238725e+01 3.246373e+00 6.748461e+00 1.007309e+01 1.212794e+01
## rhat[3] 2.647343e+02 1.351372e+01 2.387719e+02 2.555141e+02 2.645234e+02
## rhat[4] 1.034319e+01 2.891069e+00 5.393315e+00 8.284128e+00 1.010803e+01
## rhat[5] 7.895789e+01 8.566517e+00 6.302296e+01 7.298317e+01 7.863794e+01
## rhat[6] 2.385914e+02 9.878813e+00 2.195289e+02 2.317789e+02 2.384629e+02
## dev[1] 1.004898e+00 1.414663e+00 8.763011e-04 1.040125e-01 4.549952e-01
## dev[2] 8.655384e-01 1.224290e+00 1.000349e-03 8.649075e-02 3.943995e-01
## dev[3] 1.302654e+00 1.622315e+00 1.423556e-03 1.702855e-01 7.011324e-01
## dev[4] 8.329961e-01 1.178366e+00 9.743222e-04 8.366343e-02 3.804350e-01
## dev[5] 1.192679e+00 1.638148e+00 1.377183e-03 1.250219e-01 5.548485e-01
## dev[6] 1.413025e+00 1.419265e+00 4.871902e-03 3.434598e-01 9.987599e-01
## resdev 6.611791e+00 3.019223e+00 2.706447e+00 4.409680e+00 5.958335e+00
## deviance 4.697498e+01 3.019223e+00 4.306964e+01 4.477287e+01 4.632153e+01
## 75% 97.5% Rhat n.eff overlap0 f
## a 1.063410e-01 1.075857e-01 NA 1 0 1
## b 1.634924e-02 2.214316e-02 NA 1 0 1
## c 1.773858e-02 1.888898e-02 NA 1 0 1
## d 2.574387e-02 3.510430e-02 NA 1 0 1
## e 3.350695e-04 4.810107e-04 NA 1 0 1
## p[1] 1.063410e-01 1.075857e-01 NA 1 0 1
## p[2] 1.634924e-02 2.214316e-02 NA 1 0 1
## p[3] 1.773858e-02 1.888898e-02 NA 1 0 1
## p[4] 2.574387e-02 3.510430e-02 NA 1 0 1
## p[5] 6.208112e-04 7.083781e-04 NA 1 0 1
## p[6] 2.396131e-03 2.526763e-03 NA 1 0 1
## r[1] 1.104400e+04 1.104400e+04 NA 1 0 1
## r[2] 1.200000e+01 1.200000e+01 NA 1 0 1
## r[3] 2.520000e+02 2.520000e+02 NA 1 0 1
## r[4] 1.000000e+01 1.000000e+01 NA 1 0 1
## r[5] 7.400000e+01 7.400000e+01 NA 1 0 1
## r[6] 2.540000e+02 2.540000e+02 NA 1 0 1
## rhat[1] 1.112083e+04 1.125098e+04 NA 1 0 1
## rhat[2] 1.442003e+01 1.953027e+01 NA 1 0 1
## rhat[3] 2.736708e+02 2.914192e+02 NA 1 0 1
## rhat[4] 1.217685e+01 1.660433e+01 NA 1 0 1
## rhat[5] 8.451662e+01 9.643789e+01 NA 1 0 1
## rhat[6] 2.450931e+02 2.584550e+02 NA 1 0 1
## dev[1] 1.334511e+00 5.004900e+00 NA 1 0 1
## dev[2] 1.143152e+00 4.395280e+00 NA 1 0 1
## dev[3] 1.830508e+00 5.695681e+00 NA 1 0 1
## dev[4] 1.106618e+00 4.239322e+00 NA 1 0 1
## dev[5] 1.603688e+00 5.808709e+00 NA 1 0 1
## dev[6] 2.070227e+00 5.161423e+00 NA 1 0 1
## resdev 8.155144e+00 1.424058e+01 NA 1 0 1
## deviance 4.851834e+01 5.460377e+01 NA 1 0 1plot(results)
hivepi_6_xval
Simple version of the HIV Epi model with just 6 data points. (Yes, it has exactly the same comment as last example. But how is this one different?)
# Model
model_code <- "model
{
# SET PRIORS
a ~ dbeta( 1,2)
z ~ dbeta (1,1)
b <- z * (1-a) # sets constraint (1-a-b > 0)
c ~ dbeta (1,1)
d ~ dbeta (1,1)
e ~ dbeta (1,1)
# VECTOR p[1:6] HOLDS THE EXPECTED PROBABILITIES FOR EACH DATA POINT
p[1] <- a
p[2] <- b
p[3] <- c
p[4] <- d
p[5] <- (d*b + e*(1-a-b))/(1- a)
p[6] <- c*a + d*b + e*(1-a-b)
# LIKELIHOOD AND DIAGNOSTICS
for(i in 1:3) { r[i] ~ dbin(p[i],n[i])
rhat[i] <- p[i] * n[i]
dev[i] <- 2 * (r[i] * log(r[i]/rhat[i]) + (n[i]-r[i]) * log((n[i]-r[i])/(n[i]-rhat[i])))
}
for(i in 5:6) { r[i] ~ dbin(p[i],n[i])
rhat[i] <- p[i] * n[i]
dev[i] <- 2 * (r[i] * log(r[i]/rhat[i]) + (n[i]-r[i]) * log((n[i]-r[i])/(n[i]-rhat[i])))
}
dev[4]<-0
resdev <- sum(dev[])
#cross-validation of data point 4
r.rep ~ dbin(p[4],n[4])
p.xval <- step(r.rep - r[4]) - 0.5 * equals(r.rep,r[4])
}
" %>% textConnection
# Data (same as before)
data <- list(
r=c(11044, 12, 252, 10, 74, 254),
n=c(104577, 882, 15428, 473, 136139, 102287)
)
# Starting/initial values: None specified
parameters_to_save <- c("a", "b", "c", "d", "e", "p", "r", "rhat", "dev", "resdev", "p.xval")
results <- jags(data = data,
# inits = initial_values,
parameters.to.save = parameters_to_save,
model.file = model_code,
n.chains = 1, # How many chains?
n.adapt = 100,
n.iter = 50000,
n.burnin = 20000,
n.thin = 2)##
## Processing function input.......
##
## Done.
##
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 5
## Unobserved stochastic nodes: 6
## Total graph size: 93
##
## Initializing model
##
## Adaptive phase, 100 iterations x 1 chains
## If no progress bar appears JAGS has decided not to adapt
##
##
## Burn-in phase, 20000 iterations x 1 chains
##
##
## Sampling from joint posterior, 30000 iterations x 1 chains
##
##
## Calculating statistics.......
##
## Done.summary(results)## mean sd 2.5% 25% 50%
## a 1.057004e-01 9.434816e-04 1.038805e-01 1.050636e-01 1.056868e-01
## b 1.358809e-02 3.974768e-03 6.948151e-03 1.070619e-02 1.324155e-02
## c 1.715228e-02 8.648453e-04 1.545648e-02 1.656562e-02 1.715344e-02
## d 2.158887e-02 1.514051e-02 9.007164e-04 9.683381e-03 1.975353e-02
## e 2.912632e-04 1.766201e-04 1.365644e-05 1.397866e-04 2.818106e-04
## p[1] 1.057004e-01 9.434816e-04 1.038805e-01 1.050636e-01 1.056868e-01
## p[2] 1.358809e-02 3.974768e-03 6.948151e-03 1.070619e-02 1.324155e-02
## p[3] 1.715228e-02 8.648453e-04 1.545648e-02 1.656562e-02 1.715344e-02
## p[4] 2.158887e-02 1.514051e-02 9.007164e-04 9.683381e-03 1.975353e-02
## p[5] 5.851106e-04 6.389430e-05 4.664383e-04 5.410095e-04 5.828977e-04
## p[6] 2.336223e-03 9.625111e-05 2.151411e-03 2.271167e-03 2.334889e-03
## r[1] 1.104400e+04 0.000000e+00 1.104400e+04 1.104400e+04 1.104400e+04
## r[2] 1.200000e+01 0.000000e+00 1.200000e+01 1.200000e+01 1.200000e+01
## r[3] 2.520000e+02 0.000000e+00 2.520000e+02 2.520000e+02 2.520000e+02
## r[4] 1.000000e+01 0.000000e+00 1.000000e+01 1.000000e+01 1.000000e+01
## r[5] 7.400000e+01 0.000000e+00 7.400000e+01 7.400000e+01 7.400000e+01
## r[6] 2.540000e+02 0.000000e+00 2.540000e+02 2.540000e+02 2.540000e+02
## rhat[1] 1.105384e+04 9.866647e+01 1.086351e+04 1.098723e+04 1.105241e+04
## rhat[2] 1.198469e+01 3.505746e+00 6.128269e+00 9.442856e+00 1.167905e+01
## rhat[3] 2.646254e+02 1.334283e+01 2.384626e+02 2.555744e+02 2.646433e+02
## rhat[5] 7.965637e+01 8.698506e+00 6.350044e+01 7.365249e+01 7.935512e+01
## rhat[6] 2.389653e+02 9.845237e+00 2.200614e+02 2.323109e+02 2.388288e+02
## dev[1] 9.934995e-01 1.418756e+00 8.343463e-04 9.890653e-02 4.442804e-01
## dev[2] 1.057669e+00 1.524609e+00 1.103423e-03 1.021089e-01 4.862977e-01
## dev[3] 1.277279e+00 1.540366e+00 1.826850e-03 1.727566e-01 7.207671e-01
## dev[4] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## dev[5] 1.298287e+00 1.741848e+00 1.261177e-03 1.391038e-01 6.313327e-01
## dev[6] 1.361105e+00 1.385083e+00 4.040936e-03 3.168897e-01 9.505058e-01
## resdev 5.987839e+00 2.894255e+00 2.402793e+00 3.877467e+00 5.322853e+00
## p.xval 4.479000e-01 4.855766e-01 0.000000e+00 0.000000e+00 0.000000e+00
## deviance 4.221527e+01 2.894255e+00 3.863022e+01 4.010490e+01 4.155028e+01
## 75% 97.5% Rhat n.eff overlap0 f
## a 1.063193e-01 1.075926e-01 NA 1 0 1
## b 1.601359e-02 2.237791e-02 NA 1 0 1
## c 1.773533e-02 1.885851e-02 NA 1 0 1
## d 3.036678e-02 5.680448e-02 NA 1 0 1
## e 4.358797e-04 6.170680e-04 NA 1 0 1
## p[1] 1.063193e-01 1.075926e-01 NA 1 0 1
## p[2] 1.601359e-02 2.237791e-02 NA 1 0 1
## p[3] 1.773533e-02 1.885851e-02 NA 1 0 1
## p[4] 3.036678e-02 5.680448e-02 NA 1 0 1
## p[5] 6.267139e-04 7.150763e-04 NA 1 0 1
## p[6] 2.399721e-03 2.528530e-03 NA 1 0 1
## r[1] 1.104400e+04 1.104400e+04 NA 1 0 1
## r[2] 1.200000e+01 1.200000e+01 NA 1 0 1
## r[3] 2.520000e+02 2.520000e+02 NA 1 0 1
## r[4] 1.000000e+01 1.000000e+01 NA 1 0 1
## r[5] 7.400000e+01 7.400000e+01 NA 1 0 1
## r[6] 2.540000e+02 2.540000e+02 NA 1 0 1
## rhat[1] 1.111855e+04 1.125171e+04 NA 1 0 1
## rhat[2] 1.412399e+01 1.973732e+01 NA 1 0 1
## rhat[3] 2.736206e+02 2.909490e+02 NA 1 0 1
## rhat[5] 8.532020e+01 9.734978e+01 NA 1 0 1
## rhat[6] 2.454602e+02 2.586357e+02 NA 1 0 1
## dev[1] 1.310196e+00 5.060578e+00 NA 1 0 1
## dev[2] 1.392216e+00 5.282334e+00 NA 1 0 1
## dev[3] 1.830402e+00 5.563814e+00 NA 1 0 1
## dev[4] 0.000000e+00 0.000000e+00 NA 1 1 1
## dev[5] 1.751668e+00 6.159568e+00 NA 1 0 1
## dev[6] 1.969325e+00 4.995338e+00 NA 1 0 1
## resdev 7.366281e+00 1.337691e+01 NA 1 0 1
## p.xval 1.000000e+00 1.000000e+00 NA 1 1 1
## deviance 4.359371e+01 4.960434e+01 NA 1 0 1plot(results)
